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Abstract — The vestibulo-ocular reflex (VOR) consists of two 
intermingled non-linear subsystems; namely, nystagmus and 
saccade. Typically, nystagmus is analysed using a single suf- 
ficiently long signal or a concatenation of them. Saccade infor- 
mation is not analysed and discarded due to insufficient data 
length to provide consistent and minimum variance estimates. 

This paper presents a novel sparse matrix approach to 
system identification of the VOR. It allows for the simultaneous 
estimation of both nystagmus and saccade signals. We show via 
simulation of the VOR that our technique provides consistent 
and unbiased estimates in the presence of output additive noise. 

I. INTRODUCTION 

A bifurcating or switched system is one that switches 
between various modes of operation via internal or external 
influences. When a switch occurs from one subsystem to 
another, a discontinuity may result followed by a smooth 
progression under the new state [1], 

The vestibulo-ocular reflex (VOR) is a well known exam- 
ple of a biological system that exhibits non-linear bifurcating 
behavior [2], Although typically analysis of VOR dynamics 
is accomplished by relying on linear a priori modelling 
techniques, such models do not account for the rich dynamic 
behavior due to non-linearities [3], When studying subjects 
known to have vestibular deficit such models are of limited 
use due to mode interactions through initial conditions. 

Recently, Kukreja et al. [4] showed the VOR can be 
accurately modelled by the NARMAX (Non-linear AutoRe- 
gressive. Moving Average exogenous) structure. NARMAX 
models describe the current system output in terms of past 
inputs, outputs and errors and may include a variety of non- 
linear terms and, as such, they are suitable to describe the 
input-output relationship of non-linear switched systems [5], 
Moreover, Kukreja et al. [4] developed a modified extended 
least-squares (MELS) algorithm to estimate unbiased param- 
eter values of non-linear Hammerstein structure multimode 
systems. This approach uses multiple short data segments of 
both nystagmus and saccade data, which requires the signal 
be classified into modes of operation [6]. The classified 
signal is concatenated, with appropriate modelling of bifur- 
cation points, providing parameter estimates of each category 
of eye movement. This two-step process adds overhead 
to data analysis, which can be significant for sufficiently 
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long data records and batch processing. However, it can be 
reduced to one step by exploiting the nature and efficiency 
of sparse matrices. 

A sparse matrix or signal is characterised as one whose 
values are dominated by zeros and a few nonzero elements. 
A sparse structure offers the advantage of being easily 
compressed. By only storing non-zero elements a significant 
reduction in memory allocation can be realised [7]. This 
feature allows the implementation of special computational 
techniques to take advantage of the large number of zeros 
to speed up computations [7], [8]. In fact, some large sparse 
matrices do not lend themselves to efficient manipulation 
by standard algorithms designed for dense structures [9], 
[10], In this paper, we show the utilisation of a sparse 
matrix approach is ideally suited for system identification 
of switched systems with application to the VOR. 

Understanding the VOR is important for NASA’s Human 
Research Program to address (i) Risk of Impaired Con- 
trol of Spacecraft Associated Systems and (ii) Immediate 
Vehicle Egress Due to Vestibular/Sensorimotor Alterations 
Associated with Space Flight [11]. The ability to objectively 
quantify mechanisms that control VOR will support NASAs 
human exploration missions by providing methods to predict 
spatial disorientation and develop individualised countermea- 
sures for astronaut crews that must perform critical space 
operations under varied gravito-inertial conditions such as 
launch, landing and orbital maneuvering [12], 

The organisation of this paper is as follows. Section 
II introduces the NARMAX model structure. Section III 
introduces the VOR system, which is well known to ex- 
hibit bifurcating behaviour. A non-sparse matrix approach, 
namely, the MELS algorithm is summarised in §IV whilst 
a sparse matrix approach is developed in §V. This sparse 
matrix formulation is capable of simultaneously quantifing 
both nystagmus and saccade data. Section VI provides results 
of the proposed algorithm on a simulated VOR model. 
Section VII summarizes the conclusions of our study. 

II. NARMAX STRUCTURE 

A NARMAX structure’s input-output relationship is mod- 
elled as a non-linear difference equation of the form [13]. 

y(n) =f\y{n- l),...,y(n-n y ), (1) 

u{n — n u ),e(n — 1), . . -,e(n — n e )\ +e(n) 

where f is a non-linear mapping, u is the exogenous input, 
y is the output, and e is the innovation, or uncontrolled input. 

The derivation of the NARMAX model (Eqn. 1) was based 
on the zero-initial-state response. However, there are cases 



where it may be required to extend it to the non-zero- 
initial-state case [5]. This approach was taken to extend the 
NARMAX formulation to model bifurcating systems as [4] 


y m (n) =f , [ym(n~ 1), . . . ,y m (n - n y ),u m (n), . . . , ( 2 ) 

u m (n -n u ),, 8 m {n -n 5 ), e m {n- 1 ),..., 

e m (n — n e )\ +e m (n) for m = 1 , 2 ,...,.# V finite .M 


where m represents the modes of operation, u m , y m , e m are as 
defined previously, and 8 m are Kronecker impulse functions 
applied at the beginning of each data segment. Notice the 
lagged impulse values account for non-zero-initial-states. 
This modified description of the NARMAX structure allows 
it to model non-linear bifurcating systems. 

III. MODEL OF VOR DYNAMICS 

VOR dynamics are a result of ocular responses to head 
perturbations. These responses are classified as slow or fast, 
according to their average speed characteristics. The VOR 
has a time record response with a saw wave pattern called 
ocular nystagmus (Fig. 1) [3], [14], This saw wave pattern is 

Eye Position 


Output Mode2 



0 1 2 3 4 5 6 

Time (s) 

Fig. 1. VOR output. Blue Segments: slow-phase. Red Segments: fast-phase. 


a result of the dual mode switching nature of the VOR. The 
saw wave pattern represents the slow-phase which stabilizes 
the eye in space and the fast-phase which re-orients the 
eye in the direction of head rotation [3], [15], [16]. Fig. 2 
shows a typical output of the VOR illustrating its dual mode 
switching character. 



Fig. 2. Dual mode system describing vestibular dynamics. Switch position 
Si: slow-phase (nystagmus). Switch position S 2 : fast-phase (saccade) 


A. Continuous-Time Representation of VOR Model 
In Fig. 2 let f'(-), Ti(s) and Y 2 (s) are defined as [4] 
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where Y\ (s) and Y 2 (s) are first-order approximations for 
slow and fast-phase modes, Yu and Y 2 t represent the initial 
conditions and q,r are the number of switches. 

B. Discrete-Time Representation of VOR Model 

A NARMAX description of VOR slow-phases, y\ (n), and 
fast-phases, y 2 {n) of the model in Eqn. 3 is [4] 

y] (n) Dynamic Mode Si 
y 2 (n) Dynamic Mode S 2 
Pi+fhyi{n ~ 1) +J3 i[u(n)+u(n - 1)] 

+ p 4 [u 2 (n) + u 2 (n — 1)] + J85 [m 3 (n) m 3 (n — 1)] 

+ Ki 8 n(n-j)-\ \-KiSu{n- ji)\ i=l,---,q 

y 2 (n) = #1 +- 8 2 y 2 (n- 1) + 1?3 [«(«) + M (« ~ 1)] 

+ #4 [u 2 (n)+u 2 (n — 1)] + ^[u 3 (n) + u^ (n — 1)] 

+ X\ 8 u{n-k)-\ l-A ( 8 u{n-k()\ £=l,---,r. 

where 8 is the Kronecker impulse function, j,k are the lags 
on the t)|,th and 52 / th impulse and q. r are the number of 
data segments of sub-system one and two, respectively. 

Table I shows the relationship of these discrete-time (DT) 
parameters to the underlying continuous-time (CT) parame- 
ters in Eqn. 3. 
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TABLE I 

Discrete-time relationship of NARMAX model parameters to 

UNDERLYING CONTINUOUS-TIME PARAMETERS. 

IV. NON-SPARSE MATRIX APPROACH 
The steady-state response of a system can be estimated 
using the extended least-squares (ELS) algorithm [17]. How- 
ever, ELS cannot be applied to switched systems since it does 
not account for the transients which occur at each switch time 
and as a result will produce biased estimates. 

The MELS algorithm corrects for bias by including 
columns in the regressor matrix to account for initial condi- 
tions as 


* = [^1 Vs]- 


(5) 


<t> is a partitioned regressor matrix where is a function 
of input-output and prediction errors. The extension, 'Pg, 
represents the effects of initial conditions when a switch 
occurs. This yields a least-squares formulation for nystagmus 
and saccade response as 

Z\ = t f , i0i and Z 2 = < t > 202 (6) 

where Zj 2 is the measured output, 0 \ .2 the unknown system 
parameters. Notice the expression given in Eqn. 6 is a non- 
sparse matrix approach and, as such, both fast and slow phase 
need to be analysed separately. 

V. SPARSE MATRIX APPROACH 

To simplify parameter estimation and increase compu- 
tational efficiency a sparse matrix approach can be used. 
The expressions shown in Eqns. 4-6 readily lend themselves 
to a sparse matrix formulation [18]. Simultaneous analysis 
of both nystagmus and saccade signals is accomplished 
by constmcting a diagonal block-oriented data matrix and 
exploiting its natural sparseness as 
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or more compactly as 

Z = <P0. (8) 

A solution to the extended parameter set in Eqn. 8 is 

0 = ($ T it>y l & T Z. (9) 

This novel analysis technique achieves efficient computation 
of slow and fast-phase signals concurrently. 

VI. SIMULATIONS & RESULTS 

The accuracy of the sparse matrix parameter estimation 
algorithm was validated by simulating the VOR model (Fig. 
2) in continuous-time using Simulink. The parameters used 
in the simulation were typical values found in experiments 
and are shown in Table II [3]. 
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Ki 
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k 2 
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TABLE II 


Left: Continuous-time coefficient values, n : nystagmus 

TIME-CONSTANT, T 2 : SACCADE TIME-CONSTANT, K \ : NYSTAGMUS GAIN, 

K 2 : saccade gain and T: sampling interval. Right: Coefficient 

VALUES OF STATIC NONLINEARITY, a : DC TERM, b: LINEAR TERM, c: 
SQUARED TERM AND d CUBIC TERM. 

One thousand Monte-Carlo simulations were generated in 
which the input-output realisation was the same but had a 
unique Gaussian white, zero-mean, noise sequence added to 
the output. Starting with the noise free (NF) case the signal- 
to-noise ratio (SNR) of the noise sequence was decreased 
from 20 — 0 dB in increments of 5 dB. The system was 
perturbed using a sinusoid input (1/6 Hz frequency and 188 


deg/s amplitude). A sinusoid was used as input because it 
is the type of input used in a typical clinical setting. The 
system input-output was sampled at 600 Hz. 

A. Continuous-Time Parameter Estimation 

We estimated continuous-time parameters of slow and fast- 
phase dynamics using the theoretical relationships in Table 
III. Since it is impossible to measure the signal at the output 
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TABLE III 


DT TO CT RELATIONSHIPS FOR PARAMETERS T v , AND G v b. 

of the static nonlinearity, we consider the linear system to 
have unity gain and translate the overall gain onto the static 
nonlinearity giving an estimate as a product of the linear 
system and static nonlinearity gain. We deem that the best 
estimate of the linear system gain is the product of the linear 
system gain and linear coefficient of the static nonlinearity, 
i.e. G v b for v = 1 or2. 

The results of this study are shown in Fig. 3. Panels 3(a) 
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Fig. 3. Ordinate: STD about mean. Abscissa: Output SNR = NF, 20, 15, 10, 
5, 0 dB, where NF denotes noise free case. (Note that the abscissa is shown 
in decreasing SNR which corresponds to increasing noise amplitude.). 


and 3(b) display the standard deviation (STD) about the mean 
superimposed on the theoretically computed CT slow and 
fast-phase time constants. For slow-phase, the mean is close 
to the theoretical value and its STD includes the true value 
for SNR= 20- 10 dB, whilst for SNR= 5 - 0 dB the mean 
value is biased away from the true and the STD does not 
include its value. The mean of the fast-phase time constant 
is close to its true value and the STD includes the true value 
for all SNR levels. Panels 3(c) and 3(d) present the STD 
about the mean superimposed on the theoretically computed 
CT slow and fast-phase gain. The results demonstrate for all 


SNR levels the mean is close to its expected value and the 
STD encompasses the true value. 

B. Cross-Validation 

Next, we cross-validated the DT parameter estimates of 
VOR dynamics using a k-step- ahead predictor [19]. The 
quality of fit was assessed by computing the %QF as 

%QF = ) xlOO (10) 

\ N 1 / / 

where z denotes the predicted output. 

Figure 4 illustrates the predicted output for a typical 
realisation of simulated VOR using cross-validation data 
(data not used for estimation). The plot shows predicted 


Cross-Validaed Outputs QF s = 96.65, QF n = 96.26 % 



Fig. 4. Cross-validated eye position output for simulated data using a 
A-step-ahead predictor. Data noise corrupted at 15 dB SNR. 

nystagmus (red ) and saccade (green '- -’ ) superimposed 
on noise corrupted data (blue ’--’). Using our sparse matrix 
algorithm, the estimated parameters provide a slow and fast- 
phase %QF’s of 96.26% and 96.65%. 

C. Computation Time: Non-Sparse Versus Sparse Approach 

Lastly, we compared the computation time required by the 
non-sparse and sparse approach to estimate VOR parameters. 
For this study we computed the mean estimation time of 
nystagmus and saccade parameters over 1,000 trials. This 
was done to remove bias due to possible background process 
running at execution time of either approach. The analysis 
was performed in Matlab on a Macintosh with two 2.66 GHz 
6-Core Intel Xeon processors, 32 GB 1333 MHz DDR3 and 
Mac OS X 10.6.8. This study did not take advantage of the 
parallel computing available on the machine since different 
platforms have varying numbers of cores and thus, allowing 
for easy scaling for a user on his/her platform. Table IV 
presents the result of this study. The comparison shows that 


Computational Approach 

Time (s) 

Non-Sparse 

202.3 

Sparse 

16.01 


TABLE IV 

Mean computation time of 1.000 Monte Carlo simulations for 

NON-SPARSE VERSUS SPARSE MATRIX APPROACH. 

a sparse matrix approach for VOR parameter estimation is 
over 12 times faster than using a non-sparse approach. 


VII. CONCLUSIONS 

Simulation studies demonstrate that a sparse matrix ap- 
proach is a robust and efficient technique for the simulta- 
neous analysis of nystagmus and saccade dynamics. This 
method offers a tool for mathematical modelling of bio- 
logical processes, providing a quantitive measure to study 
human responses in high-g and low/zero-g environments. 
Moreover, this tool will have great significance for mission 
planners by providing them an improved ability to formu- 
late countermeasures and, hence, increasing mission safety. 
The sparse technique here allows greater insight into the 
functionality of various reflexes by providing quantitative 
measures of both fast and slow ocular dynamics from a 
single experimental record. Consequently, this approach may 
be useful to estimate the coefficients of complex non-linear 
bifurcating systems found in biology. 
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